Corcino et al.: Estimating the spread of COVID-19 


Estimation, control and forecast of COVID-19 disease spread 
in Central Visayas, Philippines 


Roberto B. Corcino*2*, Allan Roy B. Elnar+3, Gibson T. Maglasang3 
and Karl Patrick S. Casas3 
‘Research Institute for Computational Mathematics and Physics (RICMP)- 
Cebu Normal University, Cebu City, Philippines 
2Department of Mathematics, Cebu Normal University, Cebu City, Philippines 
3Department of Chemistry and Physics, Cebu Normal University, 
Cebu City, Philippines 
*Correspondence: corcinor@cnu.edu.ph 
https://doi.org/10.69721/TPS.J.2021.13.1.09 


ABSTRACT 


The current trend of infections from COVID-19 outbreak 
in Central Visayas (CV) is posing higher risk of continued 
spreading. When uncontrolled, swarming of infected individuals 
to hospitals puts a greater challenge to the health care systems in 
the region and may breakdown. It is imperative in this situation 
that data-driven decisions and policies are required the most. In 
response, this study provided estimates of the epidemiologically 
important parameters namely, reproduction metric (Ro) and 
infection (a), recovery (@) and mortality ()) rates, by using a 
modified Susceptible-Infected-Recovered-Dead (SIRD) model. 
This modified model incorporated control parameters, o and a, 
associated with enhanced community quarantine (ECQ) 
implementation and observance of social distancing (SD), 
respectively. For the covered months from 27 March 2020-10 
May 2020, results of the simulation estimated these parameters 
at: Ro=3.12, a=0.18, G= 0.029, and 7=0.029 with 90% confidence. 
Moreover, the reproduction metric can be effectively reduced 
with the combined effect of the control parameters at o,w > 0.5 
resulting to an effective Ro below unity. Similarly, this lowered 
the peak value of infection to 23% (or 7% of the total number of 
susceptible population) compared to when these measures are 
not observed and moved the peak time farther as well. While 
these estimates reflected the timely implementation of ECQ 
keeping its current level comparable with the country’s estimates 
and the world, reducing the reproduction metric effectively 
requires strict observance to both ECQ and SD control measures. 
Lastly, the temporal dynamics of this metric may not be 
necessarily true in any given area, and hence geographically 
induced. 


Keywords: Central Visayas, ECQ, social distancing, SIRD 
model, reproduction metric 
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INTRODUCTION 


The Philippines is among the hard-hit country in the ASEAN region 
from the COVID-19 outbreak which originated in Hubei, China. 
Approximately 12,400 cases have been recorded in the country as of 16 May 
2020. This is 10% compared to the overall cases recorded in the ASEAN 
region and 0.3% compared to the world data during this time. It also recorded 
2,561 recoveries and 817 deaths, hence there are 8,927 active cases in the 
entire country. 


In the local setting, Central Visayas (CV) region has been considered a 
high-risk area in the country with 1,876 active cases 53 recoveries and 28 
deaths or 1,957 recorded positive cases. With the increased turn out of active 
cases, the region has been recently placed in ECQ status quo and is constantly 
monitored by the regional interagency task force (RIATF). Putting the region 
at this quarantine level signals the inevitable increase of infections particularly 
in highly congested areas. The effects are, as Anderson et al. (2020) puts it, 
increased deaths and economic downturn are generally impacted with the 
viral spread. While measures of quarantine were put up as above, it is 
imperative in this situation that data—driven policies and decisions are needed 
in translating them into practical actions beneficial to people and the 
government. However, the challenge in the imposition of these measures may 
require stringent quarantine guidelines, such as social distancing measures 
and early implementation of quarantine, just as the Chinese government did 
when the COVID-19 outbroke in Wuhan City (Anderson et al. 2020; World 
Health Organization 2020). 


Relative to these measures, it is vital that by using epidemiological 
models, decisions have better lensing of approaches in the control and forecast 
of disease spread (Dimitrov and Meyers 2010; Heesterbeek et al. 2015; Lutz et 
al. 2019; Morgan 2019; Skrip and Townsend 2019). Thus, we provide 
estimates on epidemiologically important parameters and forecasts on the 
infection dynamics in the region using the widely-used Susceptible (S) — 
Infected (I) — Recovered (R) — Dead (D) (SIRD) model for Central Visayas 
(CV), Philippines. The bias of using the SIRD model, while it is classic, is that 
it is commonly employed by most studies of COVID-19 disease spread tackling 
the different aspects of control and forecast (Ahumada et al. 2020; 
Anastassopoulou et al. 2020; Canto and Avila-Vales 2020; Fanelli and Piazza 
2020; Shinde et al. 2020). In fact, Carletti et al. (2020) highlighted that 
COVID-19 outbreak belongs to the class of simple SIR model and its 
extensions while that of Canto and Avila-Vales (2020) showed that SIRD 
model have better predictive power over time. We have followed the method 
provided in Anastassopoulou et al. (2020) that provided reasonable picture of 
the disease dynamics at earlier stage. It is also reasonable to use the discretize 
SIRD model particularly when the element of time is on daily basis or equally 
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saying that most statistical data were obtained in discrete time (Hattaf et al. 
2015). In fact, there are several applications for the analysis of COVID-19 
spread using a discrete SIRD model and their variants including that of 
Bertozzi et al. (2020), Carcione et al. (2020), Giordano et al. (2020), Marinov 
and Marinova (2020), and ud Din et al. (2020), among others. 


We have adopted the method provided in Anastassopoulou et al. (2020), 
we determine (i.e. estimate) these essential epidemiological parameters with 
emphasis on the reproduction metric and the decay pattern of infections. 
Besides the availability of literatures to anchor on, simplicity and ease of 
algorithm, the choice of the method described therein fitted the scenario in 
Central Visayas. The urgency of putting up data-driven measures to control 
the viral spread necessitated to use the simplest algorithm from which the 
method can provide. On the same note, their method provided a coarse 
estimate of the reproduction metric without the need to compute other 
parameters, i.e. infection (a), recovery (f) and mortality (7) rates, however, 
valid only at the beginning of the disease spread. The latter justified the use 
of the model since the outbreak in Central Visayas was in its earlier stage 
during this time. 


For policy directions, we introduced control parameters associated to 
the implementation of enhanced community quarantine and observance to 
social distancing measures were introduced providing a variety of perspective 
to effectively combat and reduce the COVID-19 disease spread (Lewnard and 
Lo 2020; Sen-Crowe et al. 2020). To note, these vital information are 
important in making sound decisions as the ongoing public health crisis has 
severe social and economic repercussions. These estimates and the forecasts 
provided herein can be used as guides to local executives in planning non- 
pharmaceutical interventions to control widespread disease spreading. 


METHODS 


In this paper, we start by employing the coupled system of equations 
of the discrete SIRD model which reads 


a 


S(t) = S@-1)- e S(t — 1)/(t — 1), Eq. 1 
I(t) = 1@-1) + = S(@t-DMt-1)- BIG-1)- yIt-D, Eq.2 
R(t) = R(t—1) + BI(t— 1), Eq. 3 
D(t) = D(t-1)+ yi(t—1), Eq. 4 


where the number of susceptible, S(t), infected, I(t), recovered, R(t), and dead, 
D(t) as functions of discrete time points, t = 1, 2, 3...; and N refers to the 
population size which is a constant through time. 
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For the data resource needed in this study, we utilize the 54 data points 
for CV with reported cases of Infected (I), Recovered (R), Dead (D) from 27 
March 2020 up to 10 May 2020. For Susceptible (S) data, we used an 
approximate value of Scvy = 5M susceptible in CV. This is approximately 
equivalent to the population in Cebu Province, hence this is also a reasonable 
estimate since almost all cases recorded in CV are in Cebu Province. 


To determine whether the COVID-19 outbreak in the country and in CV 
can spread into the population or die out, we exploited the discrete method 
provided in Anastassopoulou et al. (2020) to estimate the essential 
epidemiological parameters, Ro, a, £, and ygiven by the relation 


Ro aa ae Kq. 5 

where Ro (i.e. reproduction number) represents the average number of 
secondary cases resulting in the introduction of a single infections in a totally 
susceptible population; the parameters a, £, and yare the corresponding rate 
constants associated to infected, I(t), recovered, R(t), dead, D(t) functions, 
respectively. Qualitatively, small values of Ro<1 translate to less turn out of 
secondary cases or equally saying that the spread dies out. 


It is pointed, however, that one can encounter problems of estimations 
particularly in large-scale epidemics where the actual total number of infected, 
I(t), population often is unknown. Nonetheless, Anastassopoulou et al. (2020) 
provided a coarse estimation procedure of which the researcher had faithfully 
verified. The coarse estimation of these parameters is then solved by linear 
regression using least-square problem. The parameter estimation is a two- 
step process of which rate constants @ and y are obtained using linear 
regression while a@ is obtained from the SIRD simulator. Moreover, these 
parameters are obtained for t = {1, 2, ...} with a rolling window of 1 day using 
the available daily data of actual cases from 06 March 2020 to 10 May 2020. 
Estimates and time series plots of daily Ro were rendered using Matlab 
R2019b. 


Parameter Estimation from the Reported Confirmed Cases 
We denote AI/(t), A R(t), A(t) as the reported new cases of infected, 


recovered and deaths at time ¢ respectively. The cumulative numbers of 
confirmed cases are then calculated as 


t 
CAX(t) = AX(t) 
isi 


where, X = 1,R,D, and AX = X(t) — AX(t — 1). 
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We then denote CAX(t) which is the t x 1 column vector containing 
the cumulative numbers up to time t expressed as 


CAX(t) = [CAX(1), CAX(2), ..., CAX(t)]" 


From these definitions, we can provide a coarse estimation of Ry by 
using equations 2, 3, 4, and 5. We then substitute the terms 6/(t — 1) and 
yl(t — 1) with AR = R(t) — R(t — 1) and AD = D(t) — D(t — 1) and then add 
Eqs. 4 and 5. After careful manipulation of the equations, we then arrive at 


CAI) + CAR() + CAD) _ ae _ 
CAR(t) + CAD(t) pty ”* 


Note that we can directly compute Rp with the use of regression 
without the need first to compute the other parameters. Using this equation, 
we can extract a coarse estimation of Ry by solving the linear regression using 
the least-squares problem as 


Ry = ([CAR(t) + CAD(t)]7[CAR(t) + CAD(t)])~* 
[CAR(t) + CAD(t)]"[CAI(t) + CAR(t) + CAD(t)], 


Also, the coarse estimation of the mortality and recovery rate can be 
calculated using linear regression problem for corresponding cumulative 
functions by least square given as follows 


7 = |(CAI(t) — CAD(¢) — CAR(t))' (CAI(t) — CAD(é) cAR())| 
(CAI(t) — CAD(t) — CAR(t))’ CAD(C), 


m -1 
p= [(car(e) — CAD(t) — CAR(t))' (CAI(t) — CAD(t) — CAR(t))| 
(CAI(t) — CAD(t) — CAR(t)) CAR(C). 
Next is we estimate the infection rate, a, by running the SIRD 
simulator using Band 7 as the initial conditions with one infected person 
from 27" of March until the last date available (10 May 2020). The value of a 


is then determined by wrapping around the SIRD simulator an optimization 
to solve the problem 


M 
argmin IY ewe BP)’ + woge(as BP) + wahe(aiB,9) Ip 
t=1 
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where 
fe = CAIS!R® (t) — CAI(t) 
9t = CAR*!®? (t) — CAR(t) 
h, = CADS!®°(t) — CAD(t) 


and CAX*!R? (t), (X = 1,R,D) are the cumulative cases resulting from the SIRD 
simulator at time t. Whereas the constants w,,w2, and w3 correspond to 
scalars serving in the general case as weights to the relevant functions. We can 
easily get the solution of the above problem by using the function “Isqnonlin” 
of matlab using the Levenberg-Marquard algorithm. 


SIRD Implementation 


At the onset of the spread, we set S(t=0)=N—Io, Io =1, R(t=0)=0, D(t=0) 
= 0. For the purpose of this study, we have used the percentage of the 
population instead. We implemented the model in five (5) different scenarios 
by changing the percentage of susceptible population, S(t=0), as well as 
introducing a control parameter, w, which we have associated to social 
distancing. For example, S(t=0)=N, the entire population is considered 
susceptible. With the implementation of the community quarantine, we 
assume that there is a reduction of the susceptible population using the 
relation 


n=(1-o0)N, Eq. 6 


where o=[0,1] is a constant factor related to the implementation of community 
quarantine. Qualitatively, a o=o0 means absence of ECQ implementation, 
O<o<1 means relative implementation of ECQ, and o =1 strict observance of 
ECQ. 


The last implementation of the model proposes a control parameter, 
w = [0,1], on the infection rate, a, similar to that implemented in Lin et al. 
(2010). Qualitatively, we associate this control parameter to the adherence of 
social distancing, that is an m=1 means strict social distancing is observed, 
=O means the other extreme — no observance of social distancing, and 0 < 
w < 1 means relative observance of social distancing. In control theory, the 
parameter w modulates the interaction rate of the susceptible and infected 
variables in the model, hence Eq. (2) becomes 


i(t) = it -1I) +Q-w)as(t-1i(t-1)— Bi(t-1)— yi(t—1) Eq.7 


In principle, this control parameter has also an effect on the susceptible 
function of Eq. (1), that is 


s(t) = s(t -1) -(1-o)a s(t — 1)i(@t —- 1) Eq. 8 
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with the scaled susceptible and infected function are defined as s = S/n and 
i =I/n, respectively. This is also true for Eqs. (3) and (4). We, then, make 
plots of Eq. (7) using different values of w. We note, further, that the fraction 
of susceptible (S) to the population, N, is used for the simulation. In reference 
to van den Driessche (2017), the reproduction metric in Eq. (5) becomes the 
effective reproduction number given as 


(1-o0)(1-w)a 
Ro = Eq. 
0 Bty q.9 


RESULTS 
Cumulative COVID-19 Cases in Central Visayas 


The Department of Health (DOH) updates the country of the status 
COVID-19 infections in the country. From these datasets, we used data from 
27 March 2020 to 10 May 2020 for Central Visayas data (Figure 1). It is self- 
explanatory that there is an exponential increase in the turn-out of confirmed 
cases while a linear trend is observed for recoveries and deaths. It is observed 
that during the intervals between 09 April 2020 (e.g. 36 cases) to 23 April 
2020 (e.g. 363 cases), there is a 10-fold increase in the total number of 
infections with an average of 47 daily recorded infections over a span of five 
days. 


ire data on Infected, Recovery and Death in Central Visayas 
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Figure 1. Time series data on recorded infection, recoveries, and death in 
Central Visayas, Philippines. 
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Essential Epidemiological Parameters 


From the datasets, we determine the coarse estimates of the essential 
epidemiological parameters a, /, y, and Ro. The time series of these estimates 
were plotted with a rolling window of 1 day and a 90% confidence interval. 
Using the values of these parameters, we obtain the effective reproduction 
number R, from Eq. (9) as shown in Figure 2 below. The computed Ry from 
the actual data resulted to approximately Ry ~ 3.1. It is also worth noting 
that this metric is kept at this level and is relatively comparable to the global 
metric of Ro = [2.0, 4.5]. Nonetheless, this is above the threshold value of 1. 
For purposes of comparison, the trend of the reproduction metric for the 


entire country is relatively the same, however, showing a down trend as shown 
in Figure 2d. 
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Figure 2. a) Estimated daily recovery and mortality rates (Central Visayas); b) 
estimated daily infection rate (Central Visayas); c) daily Ro values for Central 
Visayas with rolling window of 1 day from 27 March 2020 to 15 June 2020; 
and (d) daily Ro values for the Philippines (for comparison purposes) with a 
rolling window of 1 day from 11 March 2020 up to 17 April 2020. These 
estimates are obtained using linear regression at 90% confidence. 


Modified SIRD Model for Central Visayas 


The forecast on the percentage of susceptible (s), infected (i), recovered 
(r), and dead (d) individuals is shown in Figure 3. From the model, the critical 
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points we look into are peak values of these parameters and when it will occur 
when the control parameters, o and a, are implemented particularly the 
infection function, i(t). We have implemented this in different scenarios 
depicting the control of susceptible population and the observance of social 
distancing as presented below. 


It can be observed from the different scenarios in both peak values of 
infections and peak times when social distancing is introduced in the model. 
The peak values had reduced to around 23% from scenario III and V (Figure 
3). These amounts to reducing the total infections to 7% (i.e. from 1.5M to 
around 105K infections) when social distancing is observed. We also point out 
that when both ECQ and SD are strictly observed, i.e. scenario I, infections are 
reduced to almost nil and peak value occurs in a much-delayed time. We point 
out, however, that the 1.5M infected individuals is a large number to reach by 
September 2020. This large value may be caused by the computational process, 
particularly, on the assumed initial number of susceptible individuals, which 
is 5M. It would be possible that only 10% or 1% of the total population is 
susceptible to the disease because the information on how to stop the spread 
of the coronavirus has been widely disseminated throughout the country. The 
latest report in our timeline, i.e. 30 June 2020, suggested that between 33.52% 
to 40.97% of the population were observed to be susceptible corresponding to 
the number of infections of 8,660 at 90% confidence interval. 


Percentage of Total Infections forecast 
using SIRD model in Central Visayas 
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Figure 3. The infection curve simulated using modified SIRD model for 
Central Visayas in different scenarios. Thinner lines are infection curves of 
the earlier stages of the outbreak in Central Visayas. Thicker lines are 
infection curves with the inclusion of data until 15 June 2020. 
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Considering a totally susceptible population, i.e. s(t=0)=100%, 
i(t=0)=1/N, we obtained the effective reproduction metric from Eq. (9) for the 
different scenarios (Table 1). Clearly, when both ECQ and SD are observed, a 
flattened curve of infection is realized. Stricter implementations of these 
measures significantly reduce the number of infections as described by 
scenario I. 


Table 1. Qualitative description in using control metrics. 


Sample metrics used in the Modified SIRD Model 
Effective 
aeeanS euie Gea Repreducnon Description 
Metric (Ro) P 
I o=0.7 | @=0.6 Ro = 0.37 high ECQ and SD; flattened curve 
= _ a moderately high ECQ and SD; 
i BD. || Oe Ri O78 flattened the curve 
o 5 = moderately high ECQ and absence of 
ne TP pel ae Ro = 1.55 SD; approaching a flattened curve 
= Z s absence of ECQ and moderately high 
y asa eisai Ro = 1.55 SD; approaching a flattened curve 
V a 5 eas absence of ECQ and SD; not 
SO aie 0=3- flattened the curve 
DISCUSSIONS 


There are several reasons in the increasing trend of infections in CV, one 
of which is the increased turn out rate as caused by a) increased number of 
accredited testing centers, b) availability of testing kits, and c) the country- 
wide deployment of testing kits resulting to expanded testing activities. In 
perspective, these suggest that there could be more unrecorded and/or under 
reported positive cases in the country. In most pandemic studies, under- 
reporting seemed to be more likely as people experiencing mild symptoms 
often do not present themselves to health care as in the pandemic influenza 
(H1N1) (Mishra et al. 2010; Oberle et al. 2017) and TB epidemic (Zhou et al. 
2019). Mathematically, it is non-normal to observe a sudden increase/spike 
in the number of cases. With COVID-19 infections, the actual number of cases 
may be two — or three — times more since infected persons are often 
asymptomatic. This increasing trend is, in fact, not unique to Central Visayas 
and the country as a whole, rather this is true worldwide. The severity of the 
spread reflects the imposed measures to curb the outbreak in the region. To 
note, the dynamics of the spread of any disease may not necessarily be the 
same in any given area. This reflects, further, the increasing trend of the 
reproduction metric in the region compared to the reproduction number in 
the entire country, although, the estimates are approximately the same (e.g. 
Rov = 3.1 while Ro,px = 3.7). This, further, showed that the implementation 
of ECQ region-wide is relatively effective keeping the reproduction metric at 
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Ro ~ 3.10. While these values are relatively comparable worldwide and the 
earlier reported estimates in China between 2—7.1 (Lai et al. 20204; Li et al. 
2020; Mizumoto et al. 2020; Read et al. 2020; Zhang et al. 2020a; Zhou et al. 
2020), an increasing trend suggests higher risk of continued spreading. 
Consequently, this may cause social and economic instability in the region as 
the possibility of swarming of infected persons may overwhelm its health care 
facilities (Sen-Crowe et al. 2020; UPC-19PRT 2020). Moreover, our results 
complimented that of Ridenhour et al. (2018) suggesting that reproduction 
metrics varies according to geographical location affected by changes in the 
environment, population structure, viral evolution and immunity, to include 
geopolitical influences (Katz et al. 2019; Kassem 2020; Peirlinck et al. 2020). 
When these geopolitical influences cater disagreement in the implementation 
of ECQ measures, such as the earlier transition to GCQ experienced in Central 
Visayas, resulted to sudden spike in the number of infections. 


The uncertainties, however, of the actual number of infections (Lai et al. 
20204; Li et al. 2020; Zhang et al. 2020a) and of the nature and etiology of 
the virus (Andersen et al. 2020; Cascella et al. 2020; Zheng 2020) arrested the 
government and local government (LGU) leaders to implement drastic 
quarantine policies (i.e. home quarantine and social distancing) to contain the 
virus and control of the increasing loss of lives and of the continuing economic 
stress (Lim 2020). As reported in Pueyo (2020), these measures can isolate 
population to be exposed to the virus that slows down disease spreading. On 
preventing a wider spread of infection, we have shown here the effectiveness 
of how faithful adherence to the ECQ guidelines particularly on the number of 
exposed individuals as well as stricter observance to social distancing can 
reduce the reproduction metric below the threshold value. In Shim et al. 
(2020), they put prime on social distancing as measures to rapidly control the 
outbreak. This agreed well in our simulations of SD control parameter, a, 
which significantly reduced the number of infections to 7% of the initially 
assumed total susceptible population. The same observations had been 
reported in many studies of historical pandemics such as in Caley et al. (2008). 
On the contrary, Reluga (2010) pointed out that social distancing is most 
beneficial to individuals when the reproduction metric is around two and 
larger values of metric requires a more efficient social distancing measures. 
These were implied in our simulation of SD parameters in scenarios I and II 
which the effective reproduction numbers were less than unity. More 
importantly, reducing the reproduction metric effectively requires strict 
observance to both ECQ and SD control measures and/or their combined 
effects. 


We note that our predictions were proximate that of reported flattened 
curved for the number of infections in the country (Nepomuceno 2020; 
Parrocha 2020), the total number of predicted number of infections was 
relatively high given the initial number or susceptible individuals. This can be 
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justified since latest analyses suggested that it may be 10 to 40% lower where 
majority of cases were never confirmed (Annastassopoulou et al. 2020; Killeen 
and Kiware 2020; Mizumoto et al. 2020) and whose infectious carriers often 
exhibit mild to no symptoms during this time (Dong et al. 2020; Du et al. 2020; 
Lai et al. 2020b; Li et al. 2020; Tang et al. 2020; Zhang et al. 2020b; Zhou et 
al. 2020). In fact, as of the timeline of this study, the estimated susceptible 
population is between 33.52% to 40.97% corresponding to the total infections 
of 8660 at 90% confidence interval as of June 30, 2020. The latter percentage 
of susceptible population can then be used as initial input to s(t=0). Relative 
to all of these, sudden spikes in the number of infections during this time may 
be associated to contamination in densely populated or highly congested areas 
such as slum areas and/or jails (Dave et al 2020; ICRC 2020; Franco-Paredes 
et al. 2020; Simpson and Butler 2020; Wang et al 2020; WHO 2020). We 
point out, however, that most reported outbreaks in Central Visayas were 
associated to infected individuals showing mild to no symptoms at all. In 
result, granular lock down was imposed in the highly congested areas of Cebu 
City as well as several extensions of ECQ (UPCCEI 2020). 


To further improve the predictive power of the model, it is worth noting 
that granular lock down was to further control the viral spread (Basu et al. 
2020; Mandal et al. 2020; Mishra et al. 2020; Santamaria et al. 2020; 
Srivastava 2020), and hence can be a good parameter to understand the 
dynamics of COVID-19 spreading. Possible inclusion of other parameters 
such as that of density—driven dynamics (Cardoso and Goncalves 2020; 
Munshi et al. 2020), exposure—weighted dynamics (e.g. subdivision of the 
parameter I(t) into asymptomatic, mild, severe, or critical) (Adhikari et al. 
2020; Allen et al. 2020; Weitz et al. 2020) and/or age-dependent dynamics 
(Allen et al. 2020; Zhang et al. 2020b) may be seen to be a plausible direction 
to better understand the COVID-19 disease spread dynamics in Central 
Visayas or the country as a whole. 
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